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Abstract 

In  this  paper  we  study  the  superconvergence  of  the  discontinuous  Galerkin  solutions  for  nonlinear  hyperbolic  partial 
differential  equations.  On  the  first  inflow  element  we  prove  that  the  p-degree  discontinuous  finite  element  solution  con¬ 
verges  at  Radau  points  with  an  O (lf+2)  rate.  We  further  show  that  the  solution  flux  converges  on  average  at  0(h2p+2)  on 
element  outflow  boundary  when  no  reaction  terms  are  present.  For  reaction-convection  problems  we  establish  an 
0(f,m,n[2p+2'p+4))  superconvergence  rate  of  the  flux  on  element  outflow  boundary.  Globally,  we  prove  that  the  flux  con¬ 
verges  at  O (h2p+i)  on  average  at  the  outflow  of  smooth-solution  regions  for  nonlinear  conservation  laws.  Numerical 
computations  indicate  that  our  results  extend  to  nonrectangular  meshes  and  nonuniform  polynomial  degrees.  We  fur¬ 
ther  include  a  numerical  example  which  shows  that  discontinuous  solutions  are  superconvergent  to  the  unique  entropy 
solution  away  from  shock  discontinuities. 
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1.  Introduction 

The  discontinuous  Galerkin  (DG)  finite  element  method  has  been  used  to  solve  first-order  hyperbolic 
problems  and  is  gaining  in  popularity.  The  DG  method  was  first  used  for  the  neutron  equation  [24].  Since 
then,  DG  methods  have  been  used  to  solve  hyperbolic  [7-10,15,14,16,20],  parabolic  [17,18],  and  elliptic 
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[6,5,25]  partial  differential  equations.  For  a  more  complete  list  of  citations  on  the  DG  methods  and  its 
applications  consult  [13].  A  main  advantage  of  using  discontinuous  finite  element  basis  is  to  simplify  adap¬ 
tive  p-  and  /z-refinement  with  hanging  nodes. 

High-order,  p  >  0,  DG  solutions  for  nonlinear  hyperbolic  problems  exhibit  spurious  oscillations  near 
discontinuities.  These  oscillations  may  be  reduced  by  using  either  limiting  [10,11]  or  shock  capturing 
[12,21]  techniques  that  force  the  DG  solution  to  converge  to  the  unique  entropy  solution  under  mesh  refine¬ 
ment.  Many  techniques  to  suppress  spurious  oscillations  have  been  suggested  but  none  is  totally  successful. 
From  previous  computational  experience  [1],  we  discovered  that  limiting  can  reduce  spurious  oscillations 
near  shock  discontinuities  but  does  not  enhance  superconvergence  properties  of  the  DG  solution  near 
shocks.  For  these  reasons,  we  restrict  our  superconvergence  error  analysis  to  the  local  error  behavior  in 
smooth-solution  regions. 

Recently,  Adjerid  et  al.  [1]  proved  that  smooth  DG  solutions  of  one-dimensional  linear  and  nonlinear 
hyperbolic  problems  using  p-degree  polynomial  approximations  exhibit  an  O (hp+2)  superconvergence  rate 
at  the  roots  of  Radau  polynomial  of  degree  p+  1 .  They  used  this  result  to  construct  asymptotically  correct 
a  posteriori  error  estimates.  They  further  established  a  strong  0(h2p+i)  superconvergence  at  the  downwind 
end  of  every  element.  Krivodonova  and  Flaherty  [22]  proved  a  superconvergence  result  on  average  on  the 
outflow  edge  of  every  element  of  unstructured  triangular  meshes  and  constructed  a  posteriori  error  esti¬ 
mates  that  converge  to  the  true  error  under  mesh  refinement.  Adjerid  and  Massey  [4]  extended  these  results 
for  multi-dimensional  problems  using  rectangular  meshes  and  presented  an  error  analysis  for  linear  prob¬ 
lems  and  problems  with  a  nonlinear  reaction  term.  They  showed  that  the  leading  term  in  the  true  local  error 
is  spanned  by  two  [p  +  l)-degree  Radau  polynomials  in  the  x-  and  ^-directions,  respectively.  They  further 
showed  that  a  p-degree  discontinuous  finite  element  solution  exhibits  an  0(W,+2)  superconvergence  at 
Radau  points  obtained  as  a  tensor  product  of  the  roots  of  Radau  polynomial  of  degree  p  +  1.  For  a  linear 
model  problem  they  established  that,  locally,  the  solution  flux  is  O (h2p+2)  superconvergent  on  average  on 
the  outflow  element  boundary  and  the  global  solution  flux  converges  at  an  O {hLp+l)  rate  on  average  at  the 
outflow  boundary  of  the  domain.  They  used  these  superconvergence  results  to  construct  asymptotically 
exact  a  posteriori  error  estimates  for  linear  and  nonlinear  hyperbolic  problems.  In  this  paper,  we  extend 
the  error  analysis  of  Adjerid  and  Massey  [4]  to  nonlinear  hyperbolic  scalar  problems  of  the  form 


V-F (u)  =  h(x,y),  (x,y)  <E  0  =  [0,  l]2 

(11) 

and 

V  ■  F(«)  +  <f>(u)  =  h(x,y),  (x,y)  €  0  =  [0,  l]2, 

(1.2) 

with  boundary  conditions 

«Lin  =  s- 

(1.3) 

The  inflow  and  outflow  boundaries  are  defined  as 

0Oin  =  j(x,.y)  €  00  ^  •  v  <  0  j 

(1.4a) 

and 

0Oou,  =  j  (*,.)')  e  00  ^  •  V  >  oj, 

(1.4b) 

where  the  boundary  of  0,  00  =  0Ojn  U  0Oout  and  v  is  the  outward  unit  normal  to  00.  The  difficulty  with 
nonlinear  conservation  laws  (1.1)  is  that  in  general  for  smooth  flux  function  F(x,y)  and  smooth  boundary 
conditions  g,  smooth  solutions  do  not  in  general  exist  for  all  (x,y)  e  O.  Thus,  only  weak  solutions  can  be 
defined.  Furthermore,  a  weak  solution  is  unique  if  it  satisfies  the  entropy  condition  [19]. 
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In  order  to  perform  an  error  analysis  on  the  first  inflow  element,  we  assume  F  :  R  — *  R2,  <f> :  R  — >  R, 
u  :  R2  — >  R,  h  and  g  to  be  analytic  functions.  On  the  first  inflow  element  we  show  that  the  DG  solution 
of  (1.1)  is  O (ff+2)  superconvergent  at  Radau  points  and  the  leading  term  in  the  error  is  a  linear  combina¬ 
tion  of  two  Radau  polynomials.  Moreover,  the  flux  is  O (h2p+2)  superconvergent  on  average  on  the  outflow 
boundary  of  the  first  inflow  element.  For  reaction  problems  (1.2)  the  flux  is  0(/zm,n(2y,+2j>+4))  superconver¬ 
gent  on  the  outflow  boundary  of  the  first  inflow  element. 

If  we  further  assume  that  u(x,y)  is  smooth  on  Q  such  that 

Q  C  Q  and  dQin  C  3 Qin,  (1.5) 

then  the  flux  is  O (h2p+l)  superconvergent  on  average  at  the  outflow  boundary  3 Qoat. 

Finally,  computational  results  for  a  problem  with  a  shock  discontinuity  reveal  that  similar  local  super¬ 
convergence  results  hold  in  smooth-solution  regions  away  from  the  shock.  Thus,  the  error  in  smooth-solu¬ 
tion  regions  propagates  at  a  higher  order. 

This  paper  is  organized  as  follows:  In  Section  2  we  state  and  prove  the  main  superconvergence  results.  In 
Section  3  we  show  numerical  results  for  several  problems.  We  conclude  with  a  few  remarks. 


2.  Error  analysis 


In  this  section  we  will  analyze  the  DG  discretization  error  and  show  that  the  leading  term  in  the  error  is 
proportional  to  (p  +  l)-degree  Radau  polynomials  in  the  x-  and  y-directions.  Prior  to  proving  this  result  we 
need  to  recall  a  few  preliminary  lemmas. 

The  weak  discontinuous  Galerkin  formulation  is  obtained  by  partitioning  the  domain  Q  into  N  =  nxn 
square  elements  and  starting  the  integration  with  elements  whose  inflow  boundary  is  on  the  domain  inflow 
boundary. 

In  order  to  perform  an  error  analysis  we  consider  the  first  element  A  —  [0 ,hf  where  h  =  l/n  and  the  space 
Yp  of  polynomial  functions  such  that 

^iCfpU^1/1},  p>  0,  (2.1a) 

where  is  the  space  of  polynomials  of  degree  k 


(2.1b) 


These  spaces  are  suboptimal  but  they  lead  to  a  very  simple  a  posteriori  error  estimator.  For  efficiency  rea¬ 
sons  we  consider  the  smallest  spaces  that  satisfy  (2.1) 

rp  =  {v\v  =  J2  +  ^cf+1xy+1-f).  (2.2) 

(  k= o  ;=o  i=i  J 

We  note  that  tensor  product  elements  satisfy  (2.1). 

Assuming  ^(w(0, 0))  =  [ai,a2]<,  with  a,->0,  i  =  1,2,  one  can  prove  that  for  h  small  enough  the  inflow 
boundary  of  A  is  Fjn  =  Fi  U  where  F]  =  {(x,0),  0 < x < h}  and  F4  =  {(0,y),  0  <y<h).  The  outflow 
boundary  Fout  =  F2  U  F3  with  r2  =  {( h,y ),  0  <y  <  h]  and  F3  =  {(x,h),  0  <  x  <  h}. 

The  discontinuous  Galerkin  method  for  (1.1)  consists  of  determining  U(x,y)  e  Y p  on  A  such  that 

J  v-  (F(C/-)-F(£/))Fdor  +  J  J^[V  ■¥(U)-h(x,y)]Vdxdy  =  0,  MV  e  Yp. 


(2.3) 
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The  boundary  data  U~  on  Fin  is 

TT-,  \  (ng,  if  (x,y)eru 

U  ( x,y )  =  <  2.4 

l  ng,  if  (x,y)  e  r4) 

where  nw  is  the  ^-degree  polynomial  that  interpolates  w  at  the  roots  of  (p  +  l)-degree  right  Radau 
polynomial 

R;+l($)  =  Vi(0  -4(£).  -1  ^  ^  <  1.  (2.5) 

with  Lp  being  Legendre  polynomial  of  degree  p. 

Once  we  determine  the  solution  on  the  first  element  A  we  proceed  to  the  elements  whose  inflow  bound¬ 
aries  are  either  on  the  inflow  boundary  of  Q  or  an  outflow  boundary  of  A  and  continue  this  process  until  the 
solution  is  determined  in  the  whole  domain.  On  an  element  whose  inflow  boundary  is  not  on  the  boundary 
of  Q,  U~  is  defined  as 

U~(x,y)  =  lim  U((x,y)  -t-sv),  (x,y)  e  rin.  (2.6) 

s— »0+ 

The  discontinuous  Galerkin  solution  satisfies  the  DG  orthogonality  condition  which  is  obtained  by  multi¬ 
plying  (1.1)  by  V  €  Y p,  integrating  over  the  element  A  and  applying  Green’s  formula  to  obtain 

f  v  •  F(u)V  da  —  [  [  F(«)-VFdxdy  =  f  [  h(x,y)Vdxdy.  (2.7) 

J r  J  J a  J  J  a 

Applying  Green’s  formula  to  (2.3)  yields 

/  v  ■  F(U~)Vda  +  [  v  •  F{U)V  da  —  f  [  F{U)  ■  VVdxdy  =  f  [  h{xty)Vdxdy.  (2.8) 

J  /'in  “  rout  J  J  A  J  J  A 

Subtracting  (2.7)  from  (2.8)  we  obtain  the  DG  orthogonality  condition 


[  v  ■  (F((/-)  —  F(«))Fd<r+  f  v  ■  (F(U) -F{u))Vda 
Jr  ;n  J  r0ui 

-  J  -  F(m))  •  VVdxdy  =  0,  We  Yp. 


Using  the  mapping  of  A  —  [0 ,hf  into  the  canonical  element  A  =  [—1,  l]2  defined  by  x  =  h(  1  +  £)/2  and 
y  =  h(l  +  r\)l 2  and  «(<!;,  tj)  =  u(x(£) , y(r/))  we  obtain  the  DG  orthogonality  condition  (2.9)  on  the  canonical 
element 

j  v-{F{U~)-F{u))Vda+  [  v  •  (F(U)  —  F(u))Vdd 

J  rxn  ‘'Tout 

-  J  J^(F(U)-F(u))-\7VdZdri  =  0,  W  €  Yp.  (2.10) 

In  the  remainder  of  this  paper  we  omit  the  A  unless  we  feel  it  is  needed  for  clarity. 

Now,  we  recall  the  following  two  preliminary  lemmas. 

Lemma  2.1.  If  Qk  eYk  and  a  6  R2  satisfy 

f  otvQkVda-  [  /a- VF&d^  =  0,  WeYp,k^p,  (2.11) 

J  foul  J  J  A 


0,  k  <  p. 


Proof.  See  Adjerid  and  Massey  [4],  □ 
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Lemma  2.2.  Let  w  6  C°°( 0,  h)  and  nw  be  a  p-degree  polynomial  that  interpolates  w  at  Radau  points  on  [  0,  h] . 
Then  the  interpolation  error 

w(x(0)  -  nw(x(0)  =  £  Q-km\  (2.13a) 

k=p-\- 1 

where 

Q-P,M)  =  2^+1)^  -  &)({  -«.)•••(«-  y  =  cP+tR}+  i(0  (2.13b) 

and 

&:(0=<+,(0r*-p- 1(0,  *>/»+!,  (2-13c) 

vWt/z  r^(0  being  a  polynomial  of  degree  k. 


Proof.  See  Adjerid  and  Massey  [4],  □ 

Now  we  are  ready  to  state  the  main  result  for  nonlinear  conservation  laws. 

Theorem  2.3.  Let  u  and  U  be  the  solution  of  ( 1.1)  and  (2.3),  respectively.  Then  the  local  finite  element  error 
e=U-u,  (2.14) 


can  be  written  as 


«(&»/)=  E  **&(«,*), 

£=/H-l 

where 

Qp+l(i,r,)  =  ^R;+l^)  +  p2R;+l(r,). 

Furthermore,  at  the  outflow  boundary  of  the  physical  element  A 

[  v  ■  (F(k)  -  F(t/))  da  =  0{h2^2). 

Jr out 

If  the  solution  is  smooth  on  Q  =  satisfying  (1.5),  we  have  the  strong  superconvergence 

v  •  (F(m)  -  F(C/))d<r  =  0(h2p+l). 


L 

J  ©«out 


(2.15) 

(2.16) 

(2.17) 


(2.18) 


Proof.  The  proof  is  established  using  the  DG  orthogonality  condition  (2.10).  □ 

First  we  write  the  Taylor  series  of  F  about  u  to  obtain 

(2.i9) 

k= i 

Assuming  U~  to  be  the  interpolant  of  u  on  the  inflow  boundaries  described  in  (2.4)  and  using  (2.13)  we  see 
that  on  an  inflow  boundary  edge 

F(tr)  -  F(ti)  =  F(l)(n)(cr  -  u)  +  0{h2p+2).  (2.20) 

The  Maclaurin  series  of  F(1)(w)  with  respect  to  h  can  be  written  as 

CO 

F(,)(«)  =  f>[V, 

/= o 


(2,21a) 


(2.27) 
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The  0(1)  term  in  (2.26)  with  V=  1  yields 

w-a(sfc$)-*  (228) 

which  in  turn  leads  to  Qo  =  0.  We  note  that,  for  instance,  for  F(w)  =  [u2/ 2,  u]’  there  exists  Qo  ^  0  solution  of 
(2.28)  which  corresponds  to  a  nonphysical  DG  solution  with  e  =  0(1).  Here,  we  will  not  consider  such  non¬ 
physical  solutions. 

By  induction  the  O (hk),  k<p  +  1  leads  to 


/  v  •  QW&Vdo  -  j  j  Qkd>[o  '  VFd£d n  -  0,  VF  e  rp. 
''Tout  J  J  A 


(2.29) 


Applying  Lemma  2. 1  we  establish  that  Qk  —  0,  k  =  1 , 2, . . .  ,p. 

Following  the  same  line  of  reasoning  as  in  [4],  we  show  that  the  leading  term  Qp+\  can  be  split  as 

„  1  d p+\U-u).c  . 

®p+'~{p+  1)!  dhp+l  ^ 

where  Qp{Z,ri)  €  "V p  and 


—  QP+i  +  Qp, 


h=0 


(2.30a) 


QP+ 1  —  cp+ 1  ■ 


l 


tf+lu 


(0,0  )R+(Z)  +  cp+r, 


1 


8^'m 


1 2P+I(p+  1)!  0xP+> p+wr',p+l2 P+'(p+\)\  0yP+i 
Substituting  (2.30)  in  the  0(//H1)  term  of  the  series  (2.26)  leads  to 


(0,0  )*;+1(i7). 


(2.30b) 


/  d<r  +  [  <t>W.vQp+lVdo-  [  [  ^■VVQp+id^dr1+  [  ^  ■  vQpVda 

J  r-m  J  r out  J  J  a  J  r0ut 

-//< 


vvQpdZdt1  =  o,  werp. 


Using  (2.13)  and  (2.30b)  we  can  show  that 

/  ^]-vQ~p+lVda+  [  ^]-vQp+lVdo-  [  /■ 

y  rin  y  /’out  y  y^ 


®o]  ’  VK2p+1  d^d?7  =  0,  VFeiT,. 


(2.31) 


(2.32) 


Combining  (2.31)  and  (2.32)  with  Lemma  2.1  leads  to  =  0.  Using  (2.  30)  we  establish  (2.16). 
Using  (2.13)  we  can  show  that 

/  d<r  =  0,  VVzr2p_k,  k  =  p+l,...,2p. 

Jr,  „ 

Using  (2.33),  the  0(/iA’),  p  +  1  ^  A:  <  2p,  term  of  (2.26)  yields 

f  <t^-vQkvdo-f  f  <t>v]-vvQkd{dt1  =  o,  wer2p_k. 

j  r0Ui  y  j  a 

Testing  against  F  =  1  we  obtain 

[  ^o1  ■  vg*  da  =  0,  k  =  p+  l,...,2p, 

J  To ut 

which  establishes  (2.17). 


(2.33) 


(2.34) 


(2.35) 
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(2.36) 


(2.37) 


Next  we  prove  global  superconvergence  by  showing  that  on  every  element  A 
f  v-(F(tr)-F(w))chr+  f  v  •  (F(t/)  —  F(u))  dcr  =  0. 

v/Tin  JT, out 

Summing  over  all  elements  At  C  Q  we  obtain 

/  v-(F(£/-)-F(«))dff+  [_  v-(F(£7) -F(«))d<r  =  0. 

J  d^jn  J  0flout 

Using  (2.13)  leads  to  (2.18).  □ 

Next,  we  will  describe  similar  results  for  problems  of  the  form  (1.2)  where  the  DG  weak  formulation 
consists  of  determining  U(x,y )  e  Yp  on  A  such  that 

/  v  •  (F(U-)  -  F{U))Vda  +  /  /[V  •  F(U)  +  m  -  h(x,y))Vdxdy  =  0,  Werp.  (2.38) 

In  the  following  theorem  we  state  a  superconvergence  result  for  nonlinear  hyperbolic  problem  with  reaction 
terms. 

Theorem  2.4.  Let  u  and  U  be  the  solution  of(  1 .2)  and  (2.38),  respectively.  If u,  F  and  rj>  are  analytic  functions, 
then  the  local  error  estimates  (2.15)  and  (2.16)  hold.  Furthermore,  we  have  the  following  superconvergence 
results  on  the  first  inflow  element 


and 


I  v  ■  (F(U)  -  F (u))  =  o(h™n(-2p+Zp+4)) 

Jr  out 

/  v(F(£/)-F(«))dff+  f  [[4>(U)-4>(u)}dxdy  =  0(h2l,+2). 
J  r out  "  J  a 


If  the  solution  is  smooth  on  Q  satisfying  (1.5)  and  Q  =  (J^,  Ah  we  have  the  strong  superconvergence 
f  v  (F(£/)  -F(«))d<r+  f  -  0(w)]dcdy  =  0(h2p+]). 

J  0ftOul  "  «  fl 

Proof.  The  DG  orthogonality  condition  is 

f  v(F(l/-)-F(w))Fd<r+  I  v  •  (F(t/)  —  F(u))V do 

JrouX 

-  J  J (F(U)  —  F(«))  •  VV  +  [0(w)  —  (f>(U)]Vdxdy  —  0,  W  €  r„. 

On  the  canonical  element  [-1,  l]2  (2.42)  becomes 

[  v(F(tr)-F(«))Fda+  /  v-(F(I/)-F(«))Fd«r 
Jria  ^r0u\ 

-  J  J{F(U)-F{u))-VV  +  ^[4>{u)-4>{U)]Vd^dr1  =  d,  VKsf,. 

Applying  Taylor  series  to  (j)  about  u  we  have 

4>{u)  -  <p{U)  =  —a{u)e  -  y  (}>"{u),  a(u )  =  <j)'(u). 


(2.39) 


(2.40) 


(2.41) 


(2.42) 


(2.43) 


(2.44a) 
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The  Maclaurin  series  of  a(u)  with  respect  to  h  yields 

a(u)  =  2 g **&({, n ),  m, n)=\  ^ 


(2.44b) 


A=0 


Now  we  substitute  (2.44),  (2.22)  and  (2.25)  in  (2.43),  collect  terms  having  the  same  powers  of  h  to  obtain 


(f  vWoKdff—  /  /w0-VFd£d»7 
\«//'out  J  JA 

+  izhk(  [  V  •  W*Kd  a-  f  [  [W*  ■  VF  -  Zt-i  V ]  d<^d>7 

k= i  W  r out  J  J  a 

+  f"  hk(  f  v  •  W*  Vda  +  f  v\VkVd(j-  [  /  [W*  •  VV  -  Zk-iV)d^drA  =  0,  VF  e  ^ 
win  "'Coo,  J  Ja  J 


(2.45) 


where 


(2.46) 


and  v  •  W*  is  given  in  (2.27). 

Following  the  same  line  of  reasoning  as  in  Theorem  2.3  we  prove  (2.15)  and  (2.16)  for  problems  with  a 
nonlinear  reaction  term.  We  note  that  the  term  in  (2.44a)  involving  e2  is  higher  order  and  does  not 
contribute  to  our  leading  terms  and  that 


7  — 


o 

m—p— 1 _ 

QlQm-li 

/= 0 


if  m  <  p, 
if  m  >  p  +  1 . 


(2.47) 


We  prove  the  strong  superconvergence  result  (2.39)  for  nonlinear  hyperbolic  problems  with  reaction  terms 
by  setting  V  —  1  in  (2.45).  Using  (2.13),  (2.15)  and  (2.16),  to  obtain 


/  v  •  Q$Qp+\  da  —  0. 

•'  /out 

Setting  V—  1  in  the  O^),  k  >  p  +  1  in  (2.45)  leads  to 

f  v  •  Q^oQk  d<T  +  /  v-^>{,1!d<r+  [  [  Zk-\  d^df?  =  0. 
Jr-.„  7rou,  J  Ja 


(2.48) 


(2.49) 


Using  (2.22)  and  (2.47)  the  O (hp+2)  term  leads  to 


v  •  <S>Mqp+2  da  =  0. 


The  O (lf+3)  term  yields 

I  v  ^]Qp+3da  =  -f  [  Q0Qp+2 dtdr,, 

J  /’out  "  J  A 

which  is  not  necessarily  zero.  Thus,  we  establish  (2.39). 

Now,  using  (2.42)  with  V—  1  and  (2.22)  establishes  (2.40). 


(2.50) 


(2.51) 
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Letting  V=  1  in  (2.42)  and  summing  over  all  elements  At  c  Q  lead  to 

[  v  •  (F(t/~)  —  F(ti))  do-  +  /  v  •  (F(U)  -  F(u))dcr  +  [  (  [<t>{U)  -</>(«)]  dxdy  =  0.  (2.52) 

Applying  (2.13)  and  (2.22)  yields  (2.41)  which  completes  the  proof  of  Theorem  2.4.  □ 


3.  Numerical  examples 

We  will  consider  three  examples  to  validate  the  superconvergence  results  of  Section  2  for  smooth  and 
discontinuous  solutions. 


Fig.  1.  Zero-level  curves  of  the  error  for  Example  1  using  p  =  1,2,3, 4  (upper  left  to  lower  right).  Radau  points  are  shown  with  an  ‘x\ 
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Example  1.  We  consider  the  nonlinear  Burger’s  equation 

uy  +  uux=f(x,y),  (x,y)eQ,  (3.1a) 

where  Q  is  the  quadrilateral  P\P2P2P4  where  P\  =  (0,0),  P2  =  (1, 1),  P2  =  (2, 1)  and  P4  =  (-0.5,2).  We  se¬ 
lect  the  boundary  conditions  and  /  such  that  the  exact  solution  is 

u(x,y)  —  \/3  +  2 x2  +y2.  (3.1b) 

We  solve  this  problem  on  a  uniform  mesh  having  16  elements  with  p  —  1,2, 3, 4  and  plot  the  0-level 
curves  of  the  discontinuous  Galerkin  error  in  Fig.  1  with  Radau  points  marked  by  x. 


Fig.  2.  Zero-level  curves  of  the  error  for  Example  1  using  nonuniform  polynomial  degree  with  Radau  points  shown  with  an  1x’  (left). 
Distribution  of  the  polynomial  degree  of  the  finite  element  solution  (right). 


Fig.  3.  The  error  in  the  flux  on  the  outflow  boundary  of  the  first  element  versus  l//i  in  the  log-log  scale  for  Example  2. 
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These  results  show  that  the  solution  is  superconvergent  at  Radau  points  which  is  in  full  agreement  with 
Theorem  2.3.  Next  we  solve  (3.1)  with  nonuniform p  as  shown  in  Fig.  2.  The  results  shown  in  Fig.  2  indicate 
that  the  superconvergence  results  of  Section  2  are  still  valid  for  nonuniform  polynomial  degree  starting  with 
higher  degree  on  elements  at  the  inflow  boundary  of  the  domain  and  lower  polynomial  degree  on  elements 
at  the  outflow  boundary.  These  results  are  yet  to  be  proved  for  nonuniform  polynomial  degree. 

Example  2.  In  order  to  show  that  the  superconvergence  result  (2.39)  is  optimal,  we  consider  the  linear 
hyperbolic  problem  with  a  reaction  term 

ux  +  2uy  +  u=f(x,y),  (x,y)  €  [0,  l]2,  (3.2a) 

where  the  boundary  conditions  and  /  are  selected  such  that  the  exact  solution  is 

u(x,y)  =  (1  +x  +y)n.  (3.2b) 

We  solve  (3.2)  on  uniform  meshes  having  4,  16  and  36  square  elements  with  p  =  4  and  plot 

ViAf  [1)2]  •  v(u  —  U)da 

i  r ooi 

versus  l/h  in  Fig.  3.  As  predicted  by  Theorem  2.4,  these  results  show  an  0(/zmmf2/H  2,'0+4))  superconvergence 
rate  of  the  flux  on  the  outflow  boundary  of  the  first  inflow  element. 

Example  3.  Let  us  consider  the  homogeneous  inviscid  Burger’s  equation 

Uy  +  uux  =  0,  (x,y)  €  [-1,1]  x  [0, 1.999],  (3.3a) 

subject  to  the  initial  condition 

g0{x,  0)  =  1  +  sin(rcc)/2.  (3.3b) 

We  select  gi(0,y)  such  that  the  unique  entropy  solution  is  periodic  and  forms  a  shock  discontinuity  at  the 
point  (jj  -  1,|)  which  propagates  along  the  line  y  —  x  +  1.  We  solve  this  problem  on  meshes  having  5x5, 
14  x  10,  21  x  15,  28  x  20,  35  x  25,  42  x  30,  and  140  x  100  elements  with  p  =  1,2.  We  plot  the  zero-level 
curves  for  the  error  in  Fig.  5.  We  compute  the  maximum  errors  at  Radau  points  in  five  different  regions 


2 

1.8 
1.6 
1.4 
1.2 
>-  1 
0.8 
0.6 
0.4 
0.2 
0 

-1  -0.8  -0.6  -0.4  -0.2  0  0.2  0.4  0.6  0.8  1 

X 


Fig.  4.  Regions  1-5  for  problem  (3.3). 
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shown  in  Fig.  4  and  present  the  results  for  all  meshes  and  orders  in  Table  1.  We  also  plot  the  maximum 
errors  versus  mesh  size  for  p  =  1,2  in  Fig.  6.  We  note  that  the  convergence  rates  are  close  to  the  optimal 
O (hp+2)  superconvergence  rates  in  Regions  1,2,4  and  5,  for  p  =  1  and  p  =  2,  while  there  is  no  convergence 
in  Region  3  which  contains  the  shock  discontinuity.  These  results  are  in  agreement  with  the  superconver¬ 
gence  results  of  Theorem  (2.3)  for  the  entropy  solution  in  regions  away  from  the  shock  discontinuity.  We 
note  that  high-order,  p  >  0,  DG  solutions  are  known  to  develop  spurious  oscillations  near  shock  disconti¬ 
nuities  which  explains  the  nonconvergence  in  Region  3.  Spurious  oscillations  may  be  eliminated  using,  for 


Table  1 

Maximum  errors  at  Radau  points  for  problem  (3.3)  in  Regions  1-5 


NxM 

Region  1 

Region  2 

Region  3 

Region  4 

Region  5 

P-  1 

7x5 

2.3787e-03 

6.1577e-02 

4.4390e-01 

3.0670e-02 

4.7646e-02 

14x10 

5.7487e-05 

3.8567e-03 

4.4521e-01 

2.7064e-03 

9.8580e-03 

21  x  15 

1.7168e-05 

1.6498e-03 

4.3000e-01 

6.9151e-04 

4.4972e— 03 

28x20 

7.4312e-06 

1.0519e-03 

5.81 66e — 01 

3.1665e-04 

2.0990e-03 

35x25 

3.8739e-06 

1.2146e-04 

5.314Ie — 01 

1.7085e-04 

1 . 1 055e — 03 

42x30 

2.2684e— 06 

1.1820e-04 

5.8720e— 01 

1.0250e-04 

7.2652e— 04 

140  x 100 

6.3664e-08 

4.0439e-06 

6.0357e— 01 

3.2124e-06 

2. 1 1 23e — 05 

Rate 

2.9636 

2.4543 

-0.091848 

2.8665 

2.8549 

P  =  2 

7x5 

1.0780e-03 

3.5124e-02 

4.5072e-01 

2.3503e-02 

5.4503e-03 

14x10 

2.1984e-06 

9.6276e-03 

4.591 7e-01 

1.1295e— 03 

1 .51 27e — 03 

21x15 

2.041 6e-07 

2.0612e~03 

6.4203e-01 

4.2546e-04 

4.2145e-04 

28x20 

6.2627e-08 

1.9124e-04 

4.8639e-01 

7.7144e-05 

1.2900e-04 

35x25 

2.5654e-08 

4.461  le-05 

6.2678e-01 

4.8348e-06 

5.2180e— 05 

42x30 

1.2284e-08 

2.1228e-05 

6.1557e-01 

6.0801e-07 

2.8264e-05 

140x100 

9.7912e— 1 1 

8.1746e-09 

6.7466e-01 

5.5857e-09 

2.441  le-07 

Rate 

4.0167 

6.2070 

-0.053101 

4.8788 

3.8699 

Fig.  6.  Rates  of  convergence  for  problem  (3.3)  on  meshes  have  7x5,  14  x  10,  21  x  15,  28  x  20,  35  x  25,  42  x  30  and  140  x  100  elements 
on  Regions  1-5  with  p  =  1  (left)  and  p  =  2  (right). 
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instance,  shock  capturing  [21]  or  limiting  [10,11]  techniques.  Adjerid  et  al.  [1]  used  the  limiting  technique  of 
Biswas  et  al.  [10]  for  one-dimensional  nonlinear  problems  that  eliminates  spurious  oscillations  near  discon¬ 
tinuities  while  preserving  superconvergence  in  smooth  regions.  We  note  that  reliable  shock  detection  is  key 
to  successful  adaptive  limiting  strategies  [23], 


4.  Conclusions 

We  extended  the  results  of  Adjerid  and  Massey  [4]  to  nonlinear  conservation  laws.  We  proved  that  the 
discontinuous  Galerkin  finite  element  solution  is  0(/tp+2)  superconvergent  at  the  Radau  points.  We  also 
showed  that  locally  the  flux  is  0(h2p+2)  superconvergent  on  average  on  the  outflow  boundary  of  the  first 
inflow  element.  In  the  presence  of  reaction  terms  we  proved  similar  superconvergence  results  for  the  solu¬ 
tion  at  Radau  points  and  an  o(/imm(2/’+2,/,+4))  superconvergence  rate  for  the  flux  on  average  at  the  outflow 
boundary  of  the  first  inflow  element.  Furthermore,  we  showed  that  on  sub-domains  satisfying  (1.5)  the  sum 
of  the  flux  on  the  outflow  boundary  and  the  reaction  term  is  O (h2p+i)  superconvergent.  The  strong  super¬ 
convergence  of  the  flux  yields  superconvergence  of  the  solution  at  Radau  points  on  every  element.  As 
shown  in  Adjerid  and  Massey  [4],  these  superconvergence  results  for  discontinuous  finite  element  solutions 
may  be  used  to  construct  asymptotically  correct  a  posteriori  error  estimates  for  steering  adaptive  finite  ele¬ 
ment  methods.  Numerical  computations  of  Adjerid  and  Klauser  [3]  suggest  that  similar  superconvergence 
results  still  hold  for  local  discontinuous  Galerkin  solutions  of  convection-diffusion  problems.  The  error 
analysis  described  in  this  manuscript  has  been  extended  to  semi-discrete  DG  methods  for  transient  nonlin¬ 
ear  scalar  hyperbolic  conservation  laws  [2].  A  more  difficult  problem  is  to  extend  the  analysis  to  arbitrary 
elements  in  smooth-solution  regions  and  to  the  shock  region  with  limiting  or  shock  capturing.  We  will 
investigate  the  superconvergence  properties  of  the  shock  capturing  discontinuous  Galerkin  method 
[21,12]  for  hyperbolic  systems.  The  preliminary  work  of  Krivodonova  and  Flaherty  [22]  shows  supercon¬ 
vergence  of  the  flux  on  element  outflow  boundaries  for  general  unstructured  triangular  meshes,  however, 
no  pointwise  superconvergence  has  been  observed. 
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